function [TQ_tab] = mogi_inv_table(TQ_tab,stationlist,...
    vertchannel_station_list_amp,vertchannel_station_list_phase,...
    eastchannel_station_list_amp, eastchannel_station_list_phase, ...
    northchannel_station_list_amp, northchannel_station_list_phase,...
    output_pred_disp)

%{
locations = cell(size(stationlist));
locations{1} = struct('station','NPB', 'latlon', [19.41200 -155.28100 1115.0],'misfit_weight',1);
locations{2} = struct('station','SRM', 'latlon', [19.39531 -155.27400 1128.0],'misfit_weight',1);
locations{3} = struct('station','NPT', 'latlon', [19.41200 -155.28100 1115.0],'misfit_weight',1);
locations{4} = struct('station','OBL', 'latlon', [19.41750 -155.28400 1107.0],'misfit_weight',1);
locations{5} = struct('station','WRM', 'latlon', [19.40650 -155.30000 1163.0],'misfit_weight',1);
locations{6} = struct('station','SDH', 'latlon', [19.38950 -155.29400 1133.0],'misfit_weight',1);
locations{7} = struct('station','UWE', 'latlon', [19.42111 -155.29300 1240.0],'misfit_weight',1);
locations{8} = struct('station','UWB', 'latlon', [19.42469 -155.27800 1091.0],'misfit_weight',1);
locations{9} = struct('station','SBL', 'latlon', [19.42700 -155.26800 1084.0],'misfit_weight',1);
locations{10} = struct('station','HAT', 'latlon', [19.42300 -155.26100 1082.0],'misfit_weight',0);
locations{11} = struct('station','BYL', 'latlon', [19.41200 -155.26000 1059.0],'misfit_weight',0);
locations{12} = struct('station','KKO', 'latlon', [19.39781 -155.26600 1146.0], 'misfit_weight',1);
locations{13} = struct('station','RIMD', 'latlon', [19.39531 -155.27400 1128.0],'misfit_weight',1);
locations{14} = struct('station','PUHI', 'latlon', [19.385510 -155.251310 1079],'misfit_weight',1);
%convert lat-lon locations to UTM E-N coordinates
zone = utmzone(locations{1}.latlon(1),locations{1}.latlon(2));
utmstruct = defaultm('utm');
utmstruct.zone = zone;
utmstruct.geoid = wgs84Ellipsoid;
utmstruct = defaultm(utmstruct);
for i = 1:length(locations)
    locations{i}.loc = [0 0 0];
   [locations{i}.loc(1), locations{i}.loc(2), locations{i}.loc(3)] = mfwdtran(utmstruct,locations{i}.latlon(1),locations{i}.latlon(2),locations{i}.latlon(3)); 
end
%}

%station data
VariableNames = TQ_tab.Properties.VariableNames;
stationcell = {};
for ii = 1:length(VariableNames)
    if length(VariableNames{ii}) > 7
        if strcmp(VariableNames{ii}(end-3:end), '_HHZ')
            if strcmp(VariableNames{ii}(1:4), 'amp_')
                stationcell{end+1} = VariableNames{ii}(8:10);
            end
        end
    end
end
stations = struct('station',stationcell);

%station locations
for ii = 1:length(stations)
    if strcmp(stations(ii).station,'NPB')
        stations(ii).latlon = [19.41200 -155.28100 1115.0];
    end
    if strcmp(stations(ii).station,'SRM')
        stations(ii).latlon = [19.39531 -155.27400 1128.0];
    end
    if strcmp(stations(ii).station,'NPT')
        stations(ii).latlon = [19.41200 -155.28100 1115.0];
    end
    if strcmp(stations(ii).station,'OBL')
        stations(ii).latlon = [19.41750 -155.28400 1107.0];
    end
    if strcmp(stations(ii).station,'WRM')
        stations(ii).latlon = [19.40650 -155.30000 1163.0];
    end
    if strcmp(stations(ii).station,'SDH')
        stations(ii).latlon = [19.38950 -155.29400 1133.0];
    end
    if strcmp(stations(ii).station,'UWE')
        stations(ii).latlon = [19.42111 -155.29300 1240.0];
    end
    if strcmp(stations(ii).station,'UWB')
        stations(ii).latlon = [19.42469 -155.27800 1091.0];
    end
    if strcmp(stations(ii).station,'SBL')
        stations(ii).latlon = [19.42700 -155.26800 1084.0];
    end
    if strcmp(stations(ii).station,'HAT')
        stations(ii).latlon = [19.42300 -155.26100 1082.0];
    end
    if strcmp(stations(ii).station,'BYL')
        stations(ii).latlon = [19.41200 -155.26000 1059.0];
    end
    if strcmp(stations(ii).station,'KKO')
        stations(ii).latlon = [19.39781 -155.26600 1146.0];
    end
    if strcmp(stations(ii).station,'RIM')
        stations(ii).latlon = [19.39531 -155.27400 1128.0];
        stations(ii).station = 'RIMD';
    end
    if strcmp(stations(ii).station,'PUH')
        stations(ii).latlon = [19.385510 -155.251310 1079];
        stations(ii).station = 'PUHI';
    end
end
%convert lat-lon locations to UTM E-N coordinates
zone = utmzone(stations(ii).latlon(1),stations(ii).latlon(2));
utmstruct = defaultm('utm');
utmstruct.zone = zone;
utmstruct.geoid = wgs84Ellipsoid;
utmstruct = defaultm(utmstruct);
for ii = 1:length(stations)
   [stations(ii).locE, stations(ii).locN, stations(ii).locZ] = ...
       mfwdtran(utmstruct,stations(ii).latlon(1),stations(ii).latlon(2),stations(ii).latlon(3)); 
end

obs = zeros(3,length(stations)); %observation (station) locations
obstilt = zeros(3,4*length(stations));%observation locations for tilt
for i = 1:length(stations)
    obs(1,i) = stations(i).locE;
    obs(2,i) = stations(i).locN;
    obs(3,i) = 0;
    obstilt(:,4*(i-1)+1) = obs(:,i) + [-1;0;0];
    obstilt(:,4*(i-1)+2) = obs(:,i) + [1;0;0];
    obstilt(:,4*(i-1)+3) = obs(:,i) + [0;-1;0];
    obstilt(:,4*(i-1)+4) = obs(:,i) + [0;1;0];
end 

%lava lake location (from Chouet work)
epilatlon = [19.405402 -155.281041 100.0];
[epiloc_east, epiloc_north, ~] = mfwdtran(utmstruct,epilatlon(1),epilatlon(2),epilatlon(3));
mu = 10*10^9; %shear modulus
nu = 0.25; %poisson's ratio
%mogi location (from Kyle Anderson Geodetic inversions for 2018 eruptions)
kyle_east = 63.7 + 260740;
eastoffset = kyle_east-epiloc_east;
kyle_north = -41.95 + 2147480;
northoffset = kyle_north-epiloc_north;
depth = 2009;
centroideast = eastoffset+epiloc_east;
centroidnorth = northoffset+epiloc_north;

%greens functions
mdl = [depth centroideast centroidnorth 1]';
Ugreens = mogi_greensfunc(mdl,obs,mu,nu);
Ugreenstilt = mogi_greensfunc(mdl,obstilt,mu,nu);
G = [Ugreens(1,:).';Ugreens(2,:).';Ugreens(3,:).'];
Gtilt = [atan((Ugreenstilt(3,2:4:end)-Ugreenstilt(3,1:4:end))/2).'; atan((Ugreenstilt(3,4:4:end)-Ugreenstilt(3,3:4:end))/2).';zeros(length(stations),1)];
g = 9.807; %leave positive

%loop over events
TQ_tab.mogi_deltaV = NaN(height(TQ_tab),1);
TQ_tab.mogi_misfit = NaN(height(TQ_tab),1);
TQ_tab.mogi_mean_angle_misfit = NaN(height(TQ_tab),1);
% total_stationlist = [eastchannel_station_str(4:6)'; ...
%     northchannel_station_str(4:6)';...
%     vertchannel_station_str(4:6)'];
% total_channellist = [eastchannel_station_str(end-3:end)'; ...
%     northchannel_station_str(end-3:end)';...
%     vertchannel_station_str(end-3:end)'];
if output_pred_disp
   for ii = 1:length(vertchannel_station_list_amp)
       TQ_tab.(strcat(vertchannel_station_list_amp{ii},'_mogi')) = NaN(height(TQ_tab),1);
   end
   for ii = 1:length(vertchannel_station_list_phase)
       TQ_tab.(strcat(vertchannel_station_list_phase{ii},'_mogi')) = NaN(height(TQ_tab),1);
   end
   for ii = 1:length(eastchannel_station_list_amp)
       TQ_tab.(strcat(eastchannel_station_list_amp{ii},'_mogi')) = NaN(height(TQ_tab),1);
   end
   for ii = 1:length(eastchannel_station_list_phase)
       TQ_tab.(strcat(eastchannel_station_list_phase{ii},'_mogi')) = NaN(height(TQ_tab),1);
   end
   for ii = 1:length(northchannel_station_list_amp)
       TQ_tab.(strcat(northchannel_station_list_amp{ii},'_mogi')) = NaN(height(TQ_tab),1);
   end
   for ii = 1:length(northchannel_station_list_phase)
       TQ_tab.(strcat(northchannel_station_list_phase{ii},'_mogi')) = NaN(height(TQ_tab),1);
   end
end
for ii = 1:height(TQ_tab)
    %vector of displacements
    FU_real = [TQ_tab{ii,eastchannel_station_list_amp}' ...
        .*exp(1i*TQ_tab{ii,eastchannel_station_list_phase}');...
        TQ_tab{ii,northchannel_station_list_amp}' ...
        .*exp(1i*TQ_tab{ii,northchannel_station_list_phase}');...
        TQ_tab{ii,vertchannel_station_list_amp}' ...
        .*exp(1i*TQ_tab{ii,vertchannel_station_list_phase}')];
    FU_real = FU_real./(2*pi*1i*(1/TQ_tab.T(ii))); %convert from vel to disp
    %remove NaNs
    nanind = isnan(abs(FU_real));
    G_local = G(~nanind);
    Gtilt_local = Gtilt(~nanind);
%     stationlist_local = total_stationlist(nanind);
%     channellist_local = total_stationlist(nanind);
    FU_real = FU_real(~nanind);
    %invert for deltaV
    mogi_deltaV = (G_local + g*Gtilt_local/(2*pi*1i*(1/TQ_tab.T(ii)))^2)\FU_real;
    TQ_tab.mogi_deltaV(ii) = abs(mogi_deltaV);
    %use to get predicted disp
    FU_pred = (G_local + g*Gtilt_local/(2*pi*1i*(1/TQ_tab.T(ii)))^2)*mogi_deltaV;
    %calculate misfit as total error/total energy
    misfit = sum(abs(FU_real-FU_pred))/sum(abs(FU_real));
    TQ_tab.mogi_misfit(ii) = misfit;
%     misfit = sum(abs(sqrt((abs(FU_pred).*cos(angle(FU_pred)) - abs(FU_real).*cos(angle(FU_real))).^2 ...
%         + (-abs(FU_pred).*sin(angle(FU_pred)) + abs(FU_real).*sin(angle(FU_real))).^2))) ...
%         /sum(abs(FU_real));
%     TQ_tab.mogi_misfit(ii) = misfit;
    
    if output_pred_disp
        %add predicted vel to table
        FU_predvel = FU_pred.*(2*pi*1i*(1/TQ_tab.T(ii))); %convert from vel to disp
        FU_pred_nan = NaN(size(nanind));
        nancount = 0;
        for jj = 1:length(nanind)
            if ~nanind(jj)
                nancount = nancount+1;
                FU_pred_nan(jj) = FU_predvel(nancount);
            end
        end
        eaststart = 0;
        for jj = 1:length(eastchannel_station_list_amp)
            TQ_tab.(strcat(eastchannel_station_list_amp{jj},'_mogi'))(ii) = ...
                abs(FU_pred_nan(eaststart + jj));
            TQ_tab.(strcat(eastchannel_station_list_phase{jj},'_mogi'))(ii) = ...
                angle(FU_pred_nan(eaststart + jj));
        end
        northstart = length(eastchannel_station_list_amp);
        for jj = 1:length(northchannel_station_list_amp)
            TQ_tab.(strcat(northchannel_station_list_amp{jj},'_mogi'))(ii) = ...
                abs(FU_pred_nan(northstart + jj));
            TQ_tab.(strcat(northchannel_station_list_phase{jj},'_mogi'))(ii) = ...
                angle(FU_pred_nan(northstart + jj));
        end
        vertstart = length(eastchannel_station_list_amp) + length(northchannel_station_list_amp);
        for jj = 1:length(vertchannel_station_list_amp)
            TQ_tab.(strcat(vertchannel_station_list_amp{jj},'_mogi'))(ii) = ...
                abs(FU_pred_nan(vertstart + jj));
            TQ_tab.(strcat(vertchannel_station_list_phase{jj},'_mogi'))(ii) = ...
                angle(FU_pred_nan(vertstart + jj));
        end
    end
    
    %calculate angles
    count = 0;
    test_time = linspace(0,TQ_tab.T(ii),60);
    keepind = 1:length(nanind);
    keepind = keepind(~nanind);
    for jj = 1:length(stationlist)
        aeastind = jj;
        anorthind = length(stationlist) + jj;
        avertind = 2*length(stationlist) + jj;
        if sum(nanind([aeastind, anorthind, avertind])) == 0 %all channels valid
            count = count+1;
            eastind = find(aeastind == keepind);
            northind = find(anorthind == keepind);
            vertind = find(avertind == keepind);
            %approximate mean angle error numerically (could do algebraically later)
            % by adding up over one period
            east = abs(FU_real(eastind))*cos(1/TQ_tab.T(ii)*(2*pi)*test_time + angle(FU_real(eastind)));
            eastpred = abs(FU_pred(eastind))*cos(1/TQ_tab.T(ii)*(2*pi)*test_time + angle(FU_pred(eastind)));
            north = abs(FU_real(northind))*cos(1/TQ_tab.T(ii)*(2*pi)*test_time + angle(FU_real(northind)));
            northpred = abs(FU_pred(northind))*cos(1/TQ_tab.T(ii)*(2*pi)*test_time + angle(FU_pred(northind)));
            vert = abs(FU_real(vertind))*cos(1/TQ_tab.T(ii)*(2*pi)*test_time + angle(FU_real(vertind)));
            vertpred = abs(FU_pred(vertind))*cos(1/TQ_tab.T(ii)*(2*pi)*test_time + angle(FU_pred(vertind)));
            %calculate angles
            angles = (east.*eastpred + north.*northpred + vert.*vertpred);
            mag = (east.^2 + north.^2 + vert.^2).^0.5;
            magpred = (eastpred.^2 + northpred.^2 + vertpred.^2).^0.5;
            angles = acos(angles./(mag.*magpred));
            %calculate horizontal angles
            horangles = (east.*eastpred + north.*northpred);
            hormag = (east.^2 + north.^2).^0.5;
            hormagpred = (eastpred.^2 + northpred.^2).^0.5;
            horangles = acos(horangles./(hormag.*hormagpred));
            horangles = abs(wrapToPi(horangles));
            %use whatever phase closest to measured angle
            horangles(horangles > pi/2) = pi - horangles(horangles > pi/2);    
            if count == 1
                meanangles = mean(abs(wrapToPi(angles)));
                meanhorangles = mean(horangles);
            else
                meanangles = meanangles + mean(abs(wrapToPi(angles)));      
                meanhorangles = meanhorangles + mean(horangles);
            end
        end
    end
    if count > 0
        meanangles = meanangles/count;
        TQ_tab.mogi_mean_angle_misfit(ii) = meanangles*360/(2*pi);
        meanhorangles = meanhorangles/count;
        TQ_tab.mogi_mean_horangle_misfit(ii) = meanhorangles*360/(2*pi);
    end
end

end

